H3K9 methylation regulates heterochromatin silencing through incoherent feedforward loops

Histone H3 lysine-9 methylation (H3K9me) is a hallmark of the condensed and transcriptionally silent heterochromatin. It remains unclear how H3K9me controls transcription silencing and how cells delimit H3K9me domains to avoid silencing essential genes. Here, using Arabidopsis genetic systems that induce H3K9me2 in genes and transposons de novo, we show that H3K9me2 accumulation paradoxically also causes the deposition of the euchromatic mark H3K36me3 by a SET domain methyltransferase, ASHH3. ASHH3-induced H3K36me3 confers anti-silencing by preventing the demethylation of H3K4me1 by LDL2, which mediates transcriptional silencing downstream of H3K9me2. These results demonstrate that H3K9me2 not only facilitates but orchestrates silencing by actuating antagonistic silencing and anti-silencing pathways, providing insights into the molecular basis underlying proper partitioning of chromatin domains and the creation of metastable epigenetic variation.


INTRODUCTION
Ever since its discovery more than two decades ago, H3K9me has been regarded as the "heterochromatin mark, " which represses the transcriptional activity of underlying DNA (1,2).H3K9me is associated with heterochromatin in many eukaryotes including metazoans, fungi, ciliates, and plants and is involved in a variety of epigenetic phenomena, such as transposon silencing, imprinting, position effect variegation, and cell differentiation (3)(4)(5)(6).In some lineages of eukaryotes, such as vertebrates and plants, H3K9me is colocalized with DNA cytosine methylation (mC), which also plays a role in gene silencing (7,8).In the plant Arabidopsis thaliana, H3K9 dimethylation (H3K-9me2) and methylation at the CHG and CHH sites (H is A, T, or C; collectively referred to as mCH) are maintained through a selfreinforcing loop that is regulated by the H3K9 methyltransferases SUPPRESSOR OF VARIEGATION 3-9 HOMOLOGs (SUVH4/5/6) and DNA methyltransferases CHROMOMETHYLASEs (CMT2/3) (7).Although the mechanisms to establish and maintain these heterochromatin marks in repressed regions of the genome have been elucidated (7,9), how these marks silence transcription is largely unknown.In particular, the functions of these marks, aside from facilitating transcription silencing, have yet to be found.
Because heterochromatin formation is not a digital all-or-none phenomenon but rather a quantitative trait, different heterochromatin states have different dynamics and different genetic requirements (10)(11)(12).Furthermore, the inherently spontaneous nature of heterochromatin formation necessitates mechanisms that continuously counteract heterochromatin (10,12,13).Therefore, we focused on elucidating the mechanisms of establishing de novo heterochromatin formation and its counteraction using two experimental systems, in which we investigated the early steps of heterochromatin (i.e., H3K9me2/mCH) formation in genes and transposable elements (TEs) (14,15).Our results revealed the unanticipated role of H3K9me2/mCH in counteracting heterochromatin silencing by recruitment of the euchromatic mark H3K36me3, which we propose is an intrinsic fail-safe mechanism that prevents the uncontrolled spreading of heterochromatin.

Screening for factors that modulate transcription silencing
To elucidate the epigenome dynamics in heterochromatin silencing that is under the control of H3K9me2 and mCH in Arabidopsis, we leveraged a genetic mutant, increase in BONSAI methylation 1 (ibm1), which ectopically accumulates H3K9me2/mCH in the transcribed regions (gene bodies) of more than 3000 active genes (14,16,17).ibm1-induced H3K9me2/mCH in gene bodies causes a decrease in H3K4 monomethylation (H3K4me1) by the LSD1 family histone demethylase LDL2, which leads to transcriptional repression (Fig. 1A) (18).Accordingly, the developmental phenotypes of the ibm1 mutant are suppressed by the loss-of-function mutations of the LDL2 gene due to the recovery of H3K4me1 and the gene expression.Here, we focused on the pathway by which H3K9me2/mCH induces a loss of H3K4me1 by LDL2 in gene bodies.We noticed that the extent of the decrease in H3K4me1 in the ibm1 mutant showed great variation among the genes with hyper H3K9me2 in ibm1 (fig.S1A) and thus reasoned that there may be additional factor(s) that modulate the LDL2-mediated silencing pathway under the control of H3K9me2/ mCH ("silencer" and/or "anti-silencer"; Fig. 1A).To explore the underlying mechanisms, we searched for other chromatin features that are associated with the loss of H3K4me1 triggered by H3K9me2/ mCH.With linear regression analyses of H3K9me2/mCH-induced H3K4me1 changes and each epigenomic feature, including H3K4me, H3K36me, H3K27me, H2B ubiquitination (H2Bub), H2A/H3 variants, and RNA polymerase II (RNAPII), we found that the genes with low levels of H3K36me3 and H2Bub showed tendencies to lose H3K4me1 in response to H3K9me2/mCH (Fig. 1, B to D, and fig.S1,  B and C).On the other hand, genes with higher levels of H2A.Z and H2A.W were more prone to the loss of H3K4me1 induced by H3K9me2/mCH (Fig. 1B and fig.S1, B and C).

H3K9me2/mCH promotes H3K36me3 in gene bodies
From the above results, we speculated that H3K36me3 and/or H2Bub may protect against the loss of H3K4me1 in ibm1 (i.e., antisilencer); thus, we first analyzed H3K36me3 and H2Bub patterns in ibm1, ldl2, and ibm1 ldl2 double mutants by chromatin immunoprecipitation sequencing (ChIP-seq).Many genes that gained H3K9me2 and lost H3K4me1 also lost H3K36me3 in ibm1, and this loss of H3K36me3 was partially suppressed by ldl2 (Fig. 1, E to H, and .each dot represents a gene that accumulates h3K9me2 in ibm1.n = 3395.(D) intragenic patterns of h3K36me3 and h2Bub in Wt around genes categorized by h3K9me2 and h3K4me1 changes in ibm1.Among 3395 genes that accumulate h3K9me2 in ibm1 ("h3K9me2 up"), 722 genes showed clear decreases in h3K4me1 ("h3K4me1 down"), but others did not ("h3K4me1-stay") (18).(E) h3K36me3 [reads per kilobases per million mapped reads (RPKM)] in ibm1 and ldl2 for all genes (blue) and genes accumulating h3K9me2 in ibm1 (red).(F) intragenic patterns of h3K36me3 (top) and h3K4me1 (bottom) around genes categorized by h3K9me2 and h3K36me3 changes in ibm1.(G and H) h3K36me3 (G) and h3K4me1 (h) levels (RPKM) in Wt and mutants.each circle represents each gene in gene groups categorized by h3K9me2 and h3K36me3 changes in ibm1 [same as (F)]. in the boxplots, the centerline corresponds to the median, the notch represents the 95% confidence interval of the median, the upper and lower limits of the box correspond to the upper and lower quartiles, and the whiskers indicate the data range within 1.5× of the interquartile range (iQR).the P values are based on paired t tests.
fig.S2, A and B).This result is consistent with recent studies showing that H3K4me1 plays a role in recruiting the H3K36 methyltransferase(s) to gene bodies (19,20).Strikingly, however, a subset of genes that gained H3K9me2 also accumulated H3K36me3 (Fig. 1, E to G), while H2Bub did not show this change (fig.S2, A and B).The accumulation of the euchromatic mark H3K36me3 triggered by H3K9me2 was rather unexpected and suggests an unconventional mechanism.The hyperaccumulation of H3K36me3 was not suppressed by ldl2 (Fig. 1, E to G), suggesting that H3K9me2 and/or mCH, but not the loss of H3K4me1, triggered the deposition of H3K36me3 (Fig. 2A).
While the genes with increased H3K9me2 and decreased H3K36me3 were highly enriched in the genes transcriptionally down-regulated in ibm1, the genes with increased H3K9me2 and H3K36me3 were not (fig.S2C), suggesting a role of H3K36me3 in preventing transcriptional silencing.
To rule out the possibility that IBM1 itself excludes H3K36me3 from genes independent of the function to remove H3K9me2/mCH, we analyzed the ibm1 suvh4 double mutant that suppresses ibm1induced H3K9me2/mCH (14,17).The accumulation of H3K36me3 in ibm1 was completely suppressed by suvh4 (fig.S3), supporting the idea that H3K9me2/mCH paradoxically promotes the accumulation of H3K36me3 in gene bodies.Furthermore, the decrease of H3K36me3 in ibm1 was also suppressed by suvh4 (fig.S3).Therefore, our findings show that H3K9me2/mCH causes both accumulation and loss of H3K36me3 depending on the target genes.Dependency of H3K36me3 gain on SUVH4 strongly suggests that H3K9me2 and H3K36me3 co-occur in the same and/or surrounding nucleosomes in the same nucleus; however, how common the bivalent nucleosomes of H3K9me2/H3K36me3 are is still unknown.

ASHH3-mediated H3K36me3 protects against silencing
The genes that accumulated H3K36me3 in ibm1 tended to retain H3K4me1 (Fig. 1, F to H). Considering this result and the global anticorrelation between the H3K36me3 level and the loss of H3K4me1 in ibm1 (Fig. 1C), we hypothesized that H3K36me3 counteracts the loss of H3K4me1 mediated by LDL2 (Fig. 2A).To test this possibility and identify the H3K36 methyltransferase(s) responsible for H3K36me3 accumulation in ibm1, we made double mutants of IBM1 and four predicted H3K36 methyltransferase genes, ASHH1, ASHH2, ASHH3, and ASHH4.Strikingly, we found that ashh3 mutants strongly enhanced the developmental phenotypes of ibm1, while the other mutants did not (Fig. 2B and fig.S4A), suggesting that ASHH3 is responsible for the induction of H3K36me3 in response to H3K9me2/mCH in ibm1.
To test this hypothesis, we performed ChIP-seq for H3K4me1, H3K-9me2, and H3K36me3 in ibm1 and the ibm1 ashh3 double mutant.The developmental phenotypes and sterility in ibm1 progressively become severe across generations in response to the progressive accumulation of H3K9me2/mCH (17,21), but ibm1 ashh3 plants were already very defective in development and almost completely sterile in the first generation (Fig. 2B).Therefore, we performed ChIP-seq using the first generation of ibm1 and ibm1 ashh3.In the first generation, ~400 genes consistently accumulated H3K9me2 in both ibm1 and ibm1 ashh3 (fig.S5A).Although the number of genes accumulating H3K9me2 is much fewer in first-generation ibm1 and ibm1 ashh3, genes with increased H3K9me2 in third-generation ibm1 (Fig. 1F) showed similar change in H3K36me3 level in the first generation ibm1 compared to the wild type (WT), albeit to a lesser extent than in the third-generation ibm1 (Fig. 2C, left).However, the increase in H3K36me3 was not caused by ibm1 in the ashh3 background (Fig. 2C, right), confirming that ASHH3 is responsible for H3K36me3 in response to H3K9me2/mCH.Furthermore, H3K4me1 was decreased in ibm1 ashh3 compared to that in ibm1 in the genes that accumulate H3K9me2 and H3K36me3 in ibm1 (Fig. 2, D and E, and fig.S5, A  and B), supporting the idea that ASHH3-mediated H3K36me3 protects against the loss of H3K4me1 that is mediated by LDL2.The loss of H3K4me1 in ibm1 ashh3 (and in ibm1 to a smaller extent) was significantly correlated with the down-regulation of transcription, while the loss of H3K36me3 had a much weaker correlation (Fig. 2F).These results suggest that H3K4me1 is more critical for maintaining the transcriptional activity of genes than H3K36me3.Last, using FLAG-ASHH3 transgenic plants that complement the developmental phenotypes of ibm1 ashh3 double mutant (fig.S4B), we showed that genes that accumulate H3K9me2/mCH and H3K36me3 attract ASHH3 proteins in ibm1 (Fig. 2, G and H).We identified 497 genes that show increased ASHH3 localization in ibm1 compared to WT using two independent FLAG-ASHH3 transgenic plants (Fig. 2I), and these genes are significantly overlapped with genes that accumulate H3K-9me2 and H3K36me3 in ibm1 (Fig. 2J).Together, these results show that the accumulation of H3K9me2/mCH in gene bodies triggers ASHH3 recruitment and the consequent deposition of H3K36me3, which likely counteracts the loss of H3K4me1 mediated by LDL2 and transcriptional repression (Fig. 3A).The result that H3K36me3 was decreased in ibm1 ashh3 compared to that in ashh3 (Fig. 2, C  and E, and fig.S5, A and B) suggests that H3K4me1 also induces the deposition of H3K36me3, possibly by ASHH2, which recognizes H3K4me1 through the CW domain (19).
We also identified 2029 genes that are bound by ASHH3 in the WT background, albeit to a lesser extent than in ibm1 (Fig. 2K).These genes tend to accumulate more H3K36me3 than other genes (Fig. 2K), which corroborates our previous report showing that ASHH3 maintains H3K36me3 in the WT background (20).Because ASHH3 is not localized in TEs (fig.S6A), which have high H3K-9me2/mCH levels, and ASHH3 accumulation in ibm1 is positively correlated with H3K36me3 levels both in WT and ibm1 (fig.S6B), we speculate that ASHH3 is recruited to the chromatin with both H3K9me2/mCH and H3K36me3, where it reinforces the amount of H3K36me3.This idea is consistent with results showing that genes that originally have higher levels of H3K36me3 tend to gain H3K36me3 (and escape from silencing) in ibm1 (Fig. 1, E to G).These results suggest that genes with high levels of H3K36me3 will be protected from silencing induced by H3K9me2/mCH (discussed later).

Antagonistic LDL2-ASHH3 actions mediate the establishment of H3K9me2/mCH in TEs
The above results demonstrate that H3K9me2/mCH regulates transcription silencing through the combination of two incoherent feedforward loops (FFLs) (Fig. 3A and fig.S6C) (22).Given that the downstream components in these FFLs (H3K36me3 and H3K4me1) are modulated by methylation and demethylation mechanisms from outside FFLs (19,20,23,24), we speculate that incoherent FFLs, which H3K9me2/mCH induces, provide tunability to heterochromatin silencing and lead to different outcomes depending on the initial levels of H3K36me3 and H3K4me1; genes with originally higher H3K36me3 level attract ASHH3 and H3K36me3 and therefore protected from the loss of H3K4me1 and silencing even when H3K9me2/mCH accumulates, but genes with originally lower H3K-36me3 are more prone to be silenced because ASHH3 does not protect them.To elucidate whether this mechanism provides tunability to the silencing of TEs that are usually under the control of H3K9me2/ mCH, we monitored the establishment of silencing by using genetic crosses between the triple mutant for the H3K9 methyltransferase genes SUVH4, 5, 6 (suvh4/5/6) and the double mutant for the DNA methyltransferase genes CMT2, 3 (cmt2/3).Both suvh4/5/6 and cmt2/3 mutants lose H3K9me2 and mCH in TEs because these marks are maintained by a self-reinforcing loop (7).In their F1 hybrid, in which all the components of the H3K9me2-mCH feedback loop are present as heterozygous, H3K9me2/mCH is largely recovered (15).To examine how LDL2 regulates TE silencing, we conducted the genetic cross in backgrounds with or without the functional LDL2 gene (Fig. 3B) and performed whole-genome DNA methylation analyses in the F1 hybrids.In the ldl2 background, the establishment of mCH in F1 was attenuated in many TEs compared with that in F1 in the LDL2 background (Fig. 3, C and D, and fig.S7, A to D), suggesting that LDL2 promotes mCH establishment.Accordingly, the H3K4me1 level was higher in these TEs in F1-ldl2 than in F1-LDL2, and H3K-9me2 showed the opposite trend (Fig. 3, E and F, and fig.S7E).These results suggest that LDL2 promotes H3K9me2/mCH establishment by removing H3K4me1 and repressing transcription.This hypothesis is consistent with our previous results that showed that IBM1 removes H3K9me2 from transcribed genes and TEs [Fig.3A and ( 14)].In agreement with this idea, the transcription levels of the TEs with attenuated H3K9me2/mCH in F1-ldl2 were higher in F1-ldl2 than in F1-LDL2 (Fig. 3G).Furthermore, H3K4me3 and H3K36me3, transcriptionassociated histone marks, were also higher in F1-ldl2 than in F1-LDL2 (fig.S7, F and G).Together, our results demonstrate that LDL2 plays a role in the silencing establishment of TEs most likely by removing H3K4me1.The fact that not all TEs are attenuated in silencing in the ldl2 mutant (Fig. 3C) suggests that other partially redundant silencing pathways that are associated with H3K9me2/mCH are also involved.
We then tested the effects of ASHH3-H3K36me3, which counteracts LDL2, on mCH establishment in TEs by conducting the suvh4/5/6-cmt2/3 cross in the ashh3 background (Fig. 4A).We found 14 TEs that reproducibly had more mCH in the ashh3 background than in the ASHH3 background (Fig. 4, B to D, and figs.S8, A to D and S9), suggesting that ASHH3 counteracts the accumulation of mCH as an antisilencer in some TEs, as predicted by the model in Fig. 3A.Consistent with the elevated level of mCH in F1 in the ashh3 background, H3K-9me2 levels were also higher in these TEs in F1-ashh3 than in F1-ASHH3 (Fig. 4E and fig.S8E).These ASHH3-regulated TEs had higher H3K36me3 and H3K4me1 levels than other TEs in F1-ASHH3, and the levels of these marks were decreased in F1-ashh3 (Fig. 4, F and G, and fig.S8E).Furthermore, we found that the transcription of 9 of the 14 ASHH3-regulated TEs was significantly down-regulated [false discovery rate (FDR) < 0.05] in F1-ashh3 compared to F1-ASHH3 (Fig. 4H and fig.S8F).Although the ASHH3-regulated TEs are much fewer than the LDL2-regulated TEs, they are significantly overlapped (fig.S8G).
Together, our results demonstrate that the antagonistic actions of LDL2 that removes H3K4me1 and ASHH3 that deposits H3K36me3, which are under the control of H3K9me2/mCH, mediate the de novo establishment of transcription silencing via intricate feedback and feedforward mechanisms (Fig. 3A).

DISCUSSION
In this report, we demonstrated that H3K9me2/mCH both facilitates and impedes silencing, the latter of which is contrary to the widely believed canonical function of H3K9 methylation (2,6).However, these two seemingly paradoxical functions of H3K9me2 can explain the robust differentiation of actively transcribed genes and silent TEs (Fig. 5).When H3K9me2/mCH accumulates in transcribed regions, which is thought to occur constantly but transiently (16,25), genes with higher levels of H3K36me3 attract ASHH3 to methylate H3K36me3, which may repel LDL2 function and therefore protect H3K4me1 and transcription levels (Fig. 5, left).H3K9me2/mCH is removed by transcription-promoted IBM1 (14).On the other hand, genes with less H3K36me3 cannot attract ASHH3 in response to H3K9me2/mCH accumulation; thus, LDL2 functions to remove H3K4me1, which eventually leads to transcriptional silencing and further accumulation of H3K9me2/mCH (Fig. 5, right).We propose that this mechanism acts as a safeguard against the spurious silencing of essential genes.In the ashh3 background, only a small number of TEs showed increased H3K9me2/mCH (Fig. 4, B to E), indicating that ASHH3's anti-silencing function is more prominent in genes than in TEs.This differential act of ASHH3 in genes and TEs could be a basis for the differentiation of epigenetic status between genes and TEs.A recent report showed that certain genomic regions in human embryonic stem cells accumulate both H3K9me3 and H3K36me3 (26).Animals have multiple methyltransferases that catalyze H3K36me3, and the biological roles of some of these methyltransferases are still unknown (27).Hence, it is tempting to speculate that the analogous mechanism(s) of fine-tuning heterochromatin silencing may also exist in animals.Regarding the connection between H3K36me3 and silencing, it has been shown that in mammals, the establishment of DNA methylation in the bodies of transcribed genes is guided by the transcription-coupled accumulation of H3K-36me3, which recruits the DNA methyltransferase DNMT3B (28), although the function of DNA methylation in gene body is still elusive.In animals and yeast, the similar noncanonical roles of H3K9 methylation in promoting gene transcription have been reported as mechanisms to solve the dilemma that transcription needs to occur for producing small-interfering RNA (siRNA) or PIWI-interacting RNA, which are needed to maintain repressed heterochromatin status (29)(30)(31)(32)(33). Plants have RNA polymerase IV dedicated to producing siRNA from repressed repeat sequences with H3K9me2/mCH, so that the above dilemma can be circumvented (34,35).These mechanisms can be viewed as self-reinforcing feedback loops to stably maintain the repressed state of heterochromatin.Meanwhile, our study demonstrates a mechanism of incoherent FFLs, in which H3K9me2/mCHpromoted euchromatic mark H3K36me3 counteracts the establishment of silencing by inhibiting the demethylation of H3K4me1 (Fig. 5).Given that the deposition of H3K36me3 is coupled with the transcription process and H3K4me1, another transcription-promoted mark (19,20,23), we anticipate that the H3K36me3-mediated anti-silencing mechanism uncovered here helps create various metastable epigenetic statuses depending on the underlying DNA sequence, transcription potential, and intracellular and extracellular environment; thus, this mechanism could be a molecular basis for natural epigenetic variation (36,37).

Plasmid construction and transformation
pASHH3::3xFLAG-ASHH3-HA was constructed using NEBuilder HiFi DNA assembly Mix (New England Biolabs) and pPLV01 vector (41).The promoter region of the ASHH3 gene was amplified with the following primers: 5′-AATTCTAGTTGGAATGGGTTCTAAA TCAATCATCAAACGATAAG-3′ and 5′-AATCTCCATCGTGATCT TTGTAATCCATGTCGAAAGAATCTAAAAG-3′.The coding region until just before the stop codon was amplified with the following primers: 5′-ACAAAGATGATGATGATAAGCCAGCCAGCAAAA AGGTC-3′ and 5′-TCTGGAACATCGTATGGGTAGACAATCTC CCAGTCTTCTC-3′ and then with the following primers: 5′-GATT ACAAAGATCACGATGGAGATTACAAAGATCACGATATAGAT-TACAAAGATGATGATGATAAG-3′ and 5′-GATCCTTATGGAG TTGGGTTTTAAGCGTAATCTGGAACATCGTATGGGTA-3′.The amplified fragments were assembled into pPLV01 using NEBuilder.Fig. 5.A model of the mechanism that partitions epigenetic marks by counteracting LDL2 and ASHH3.When h3K9me2/mch accumulates in the transcribed region of genes or tes, genes/tes with high levels of h3K36me3 are kept in active chromatin status by the concerted functions of AShh3 and iBM1 (left).Meanwhile, genes/tes with lower levels of h3K36me3 cannot attract AShh3, and thus silencing is established through the removal of h3K4me1 by ldl2 (right).
For correlation analysis in Fig. 1 (B and C), we used 3395 proteincoding genes that accumulate H3K9me2 more in ibm1 than in WT (18).The Pearson correlation coefficient for each comparison between the explanatory variable (screening dataset described above) and the "decrease of H3K4me1 in ibm1" {(RPKM [WT] -RPKM [ibm1] )/RPKM [WT] } as the objective variable was calculated.Metaplots of chromatin features around genes categorized by H3K9me2 and H3K4me1 changes in ibm1 were drawn using the deepTools (48).

Chromatin immunoprecipitation sequencing
ChIP-seq for histone modifications was performed as described as enhanced ChIP-seq (eChIP-seq) in ( 49) with some modifications.A total of ~0.5 g of 2-week-old seedlings grown on MS plates was collected and frozen with liquid nitrogen, ground into a fine powder, and crosslinked with 1% formaldehyde in phosphate-buffered saline (PBS) buffer containing 1 mM Pefabloc SC (Merck), cOmplete proteinase inhibitor cocktail (Merck), and 0.3% Triton X-100, for 10 min in room temperature.After quenching formaldehyde by adding 200 mM glycine and incubating for 5 min, the fixed tissues were washed with ice-cold PBS once with centrifugation at 5000g for 5 min.The tissue pellet was dissolved by the low-salt ChIP buffer without Triton X-100 [50 mM Hepes-KOH (pH 7.5), 150 mM NaCl, 1 mM EDTA, 0.1% sodium deoxycholate, and 0.1% SDS] containing cOmplete proteinase inhibitor cocktail and sonicated with Picoruptor (Diagenode).After sonication, debris was removed by centrifugation at 20,000g for 10 min.The supernatant was transferred to a new tube, and Triton X-100 was added to make the final concentration of 1%.The sample was incubated with 1 μg of the following antibodies: α-H3 (ab1791, Abcam), α-H3K9me2 (MABI0317, Wako), α-H3K4me1 (ab8895, Abcam), α-H3K36me3 (ab9050, Abcam), and α-H2Bub (MM-0029, MediMabs).Antibody reaction was performed overnight at 4°C with rotation.The sample containing the antibody-chromatin complex was incubated with Dynabeads Protein G (Thermo Fisher Scientific) for H3K9me2 or Dynabeads M-280 Sheep anti-mouse immunoglobulin G (Thermo Fisher Scientific) for other antibodies for 2 hours at 4°C with rotation.The beads were washed once with low-salt ChIP buffer [50 mM Hepes-KOH (pH 7.5), 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, and 0.1% SDS] containing cOmplete proteinase inhibitor cocktail, two times with high-salt ChIP buffer [50 mM Hepes-KOH (pH 7.5), 350 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, and 0.1% SDS], once with ChIP wash buffer [10 mM tris-HCl (pH 8.0), 250 mM LiCl, 0.5% NP-40, 1 mM EDTA, and 0.1% sodium deoxycholate], and once with TE buffer.The chromatin was eluted by adding ChIP elution buffer [50 mM tris-HCl (pH 7.5), 10 mM EDTA, and 1% SDS] and incubating at 65°C for 15 min.Subsequently, 4 μl of Proteinase K (20 mg/ml; Thermo Fisher Scientific) was added to the sample and incubated at 55°C overnight.The immunoprecipitated DNA was purified using the Monarch PCR & DNA Cleanup Kit (New England Biolabs).The libraries for Illumina sequencing were constructed using the ThruPLEX DNA-Seq Kit (Clontech) and purified using SPRIselect Beads (Beckman Coulter).The sequencing was performed by the HiSeq X or NovaSeq X Plus sequencer (Illumina).Two independent biological replicates were performed for all of the ChIP-seq experiments.
ChIP-seq for 3xFLAG-ASHH3 was performed with some modifications from histone modifications eChIP-seq.A total of 0.8 g of seedlings was collected and crosslinked in 1% formaldehyde solution by applying vacuum for 15 min.Glycine was added to the solution and further incubated for 5 min with vacuum.Tissues were washed with water, dried with paper towels, frozen with liquid nitrogen, and kept in a −80°C freezer until use.The fixed tissue was ground into fine powder and proceeded as described above.The sonication was performed with Covaris S220 Focused-ultrasonicator (Covaris).Immunoprecipitation was conducted with 1 μg of α-FLAG antibody (F3165, Merck).
Reads were mapped on the Arabidopsis TAIR10 genome using Bowtie (46).The read count for each transcription unit was calculated using the coverage function of BEDTools (47), and then reads per kilobases per million mapped reads (RPKM) values were calculated.For extracting the genes with increased and decreased H3K36me3 in ibm1 (Fig. 1, F to H), RPKM values were calculated and compared between WT and ibm1 for each replicate.Genes with RPKM differences of more than three in both replicates were identified as increased or decreased H3K36me3.For statistical analyses such as in Fig. 1G, the averages of RPKM values between biological replicates were used.For extracting the genes with increased FLAG-ASHH3 in ibm1 (Fig. 2I) and high FLAG-ASHH3 in WT (ASHH3 bound; Fig. 2K), RPKM values from two biological replicates using two independent transgenic lines were used.Genes with increased FLAG-ASHH3 in ibm1 are selected by RPKM [FLAG-ASHH3 in ibm1] -RPKM [FLAG-ASHH3 in WT] > 5 in both transgenic lines.ASHH3-bound genes are selected by RPKM [FLAG-ASHH3 in WT]  -RPKM [WT control] > 2 in both TGs.Metaplots and heatmaps were drawn using the deepTools (48), and the results of one of two replicates were shown.Both replicates showed essentially the same results.

mRNA-seq
Total RNA was extracted from 2-week-old seedlings grown on an MS plate using the RNeasy Plant Mini Kit (QIAGEN) following the manufacturer's instructions.mRNA-seq libraries were constructed from 600 ng of total RNA using the KAPA mRNA HyperPrep Kit (Kapa Biosystems).Three independent biological replicates were analyzed for each genotype.The resulting libraries were 150-base pair (bp) paired-end sequenced by the HiSeq X or NovaSeq X Plus sequencer (Illumina).
Reads were mapped on the Arabidopsis TAIR10 genome using STAR aligner (50) with --outFilterType, BySJout; --alignSJoverhangMin, 8; --alignSJDBoverhangMin, 1; --clip3pNbases, 50; --quantMode GeneCounts parameters.The resulting per-gene counts were used for downstream analysis.For TE gene analysis (Figs. 3G and 4H), the strand that showed larger counts was chosen for each TE gene based on the sum of all samples and used for further analyses because some TE genes show higher transcription levels of antisense strand than sense strand.The differentially expressed gene analysis in Fig. 4H

Fig. 1 .
Fig. 1.H3K36, as a putative anti-silencer, is induced by H3K9me2/mCH.(A) hypothetical factors that modulate silencing triggered by h3K9me2/mch.(B) Screening for other chromatin features that are correlated with the loss of h3K4me1 triggered by h3K9me2/mch in ibm1.Pearson's R 2 (coefficient of determination) values for factors positively correlated with the loss of h3K4me1 (silencer candidates) and −R 2 values for factors negatively correlated with the loss of h3K4me1 (anti-silencer candidates) are shown.(C) correlations between h3K36me3 levels(20) (left) and h2Bub levels (44) (right) and the "decrease in h3K4me1 in ibm1" (see Materials and Methods) are shown as scatter plots, linear regression lines, and Pearson's R 2 .each dot represents a gene that accumulates h3K9me2 in ibm1.n = 3395.(D) intragenic patterns of h3K36me3 and h2Bub in Wt around genes categorized by h3K9me2 and h3K4me1 changes in ibm1.Among 3395 genes that accumulate h3K9me2 in ibm1 ("h3K9me2 up"), 722 genes showed clear decreases in h3K4me1 ("h3K4me1 down"), but others did not ("h3K4me1-stay")(18). (E) h3K36me3 [reads per kilobases per million mapped reads (RPKM)] in ibm1 and ldl2 for all genes (blue) and genes accumulating h3K9me2 in ibm1 (red).(F) intragenic patterns of h3K36me3 (top) and h3K4me1 (bottom) around genes categorized by h3K9me2 and h3K36me3 changes in ibm1.(G and H) h3K36me3 (G) and h3K4me1 (h) levels (RPKM) in Wt and mutants.each circle represents each gene in gene groups categorized by h3K9me2 and h3K36me3 changes in ibm1 [same as (F)]. in the boxplots, the centerline corresponds to the median, the notch represents the 95% confidence interval of the median, the upper and lower limits of the box correspond to the upper and lower quartiles, and the whiskers indicate the data range within 1.5× of the interquartile range (iQR).the P values are based on paired t tests.

Fig. 2 .
Fig. 2. ASHH3 mediates H3K9me2/mCH-triggered H3K36me3 to prevent the loss of H3K4me1.(A) Model depicting the h3K36me3 promotion by h3K9me2/mch, which prevents the h3K4me1 loss by ldl2.(B) Four-week-old Wt, ashh3, ibm1, and ibm1 ashh3 plants.(C and D) Kernel density plots showing the log 2 ratio of RPKM [h3K36me3 (c) and h3K4me1 (d)] between indicated samples.Gene categories are the same as Fig. 1 (F to h); "h3K9me2 up h3K36me3 up (164), " "h3K9me2 up h3K36me3 down (445), " and "h3K9me2 up h3K36me3 stay (2785)." (E) Browser views of h3K9me2, h3K4me1, and h3K36me3 in ibm1 and ashh3 around genes that have increased h3K9me2 and h3K36me3 in ibm1.(F) Scatter plots comparing h3K4me1 changes (left) and h3K36me3 changes (right), with mRnA changes between Wt and mutants.dots represent expressed genes with increased h3K9me2 in 1G ibm1 (352 genes).lines represent linear regression and shaded areas represent 95% confidence intervals.Spearman's correlation coefficient ρ is shown for each Wt-mutant comparison.*: 0.01 < P value < 0.05; ***P value < 0.0001.(G) FlAG-AShh3 localization in Wt and ibm1 plants around genes with increased h3K9me2 and h3K36me3 in ibm1 [shown in (e)].Wt on top is the nontransgenic negative control.(H) Average profiles of FlAG-AShh3 and h3K9me2 around genes categorized as (c) (top) and heatmaps around genes in "h3K9me2 up h3K36me3 up" category, sorted by FlAG-AShh3 level in ibm1 (bottom).(I) difference of FlAG-AShh3 signal between Wt and ibm1 in two independent transgenic lines, tG1 and tG2.Blue dots indicate all genes and red dots indicate 497 genes with more AShh3 accumulation in ibm1 than in Wt. (J) enrichment analysis of the genes with AShh3 accumulation in ibm1 (i) in genes categorized by h3K9me2 and h3K36me3 changes in ibm1.the color indicates the P value of the hypergeometric test.(K) Pattern of FlAG-AShh3 and h3K36me3 in Wt around 2029 AShh3-bound genes and other protein-coding genes (PcGs).

mCHH♀Fig. 3 .
Fig. 3. LDL2 facilitates TE silencing by removing H3K4me1.(A)A model depicting the h3K9me2/mch-promoted h3K36me3 by AShh3, which prevents the loss of h3K4me1 mediated by ldl2.transcribed sequences are targeted by iBM1 for h3K9 demethylation(14).therefore, the antagonistic actions of ldl2 and AShh3 are predicted to affect h3K9me2/mch dynamics through feedback regulation.(B) experimental design for analyzing the function of ldl2 in te silencing.(C) effects of ldl2 mutation on mchG (left) and mchh (right) establishment in sxc (x axis) and cxs (y axis).each dot represents a te-encoded gene ("te gene"; n = 3728), and red dots represent "ldl2-regulated te genes" (n = 194), which show lower mchG in sxc ldl2 and cxs ldl2 than in sxc LDL2 and cxs LDL2, respectively (see Materials and Methods).Averages of biological replicates are plotted.(D) Averaged profiles of mchG and mchh in Wt and four F1 lines around all te genes (left) and ldl2-regulated te genes (right).(E and F) h3K9me2 (e) and h3K4me1 (F) patterns in Wt and four F1 lines around all te genes (top) and ldl2-regulated te genes (bottom).(G) mRnA levels [log 2 of transcripts per million (tPM) + 0.05] of ldl2-regulated te genes in Wt, parental mutant plants, and F1 plants.each circle represents each te gene. in the boxplots, the centerline corresponds to the median, the notch represents the 95% confidence interval of the median, the upper and lower limits of the box correspond to the upper and lower quartiles, and the whiskers indicate the data range within 1.5× of the iQR. the P values are based on paired t tests.n.s., not significant.

♀Fig. 4 .
Fig. 4. ASHH3 counteracts the silencing of several TEs.(A) experimental design for analyzing the function of AShh3 in te silencing.(B) effects of ashh3 mutation on mchG (left) and mchh (right) establishment in sxc (x axis) and cxs (y axis).each dot represents a te gene (n = 3551), and red dots represent AShh3-regulated te genes (n = 14), which show higher mchG in sxc ashh3 and cxs ashh3 than in sxc ASHH3 and cxs ASHH3, respectively (see methods).Averages of biological replicates are plotted.(C and D) Averaged profiles of mchG (c) and mchh (d) in Wt and four F1 lines around all te genes (left) and AShh3-regulated te genes (right).(E to G) h3K9me2 (e), h3K36me3 (F), and h3K4me1 (G) patterns in Wt and four F1 lines around all te genes (top) and AShh3-regulated te genes (bottom).(H) volcano plot of mRnA-seq comparing tPM of cxs ashh3 and cxs ASHH3.Blue dotted lines represent FdR = 0.05.Blue dots represent all genes including protein-coding genes and te genes (n = 33,603), and red dots represent AShh3-regulated te genes.
and fig.S8F was performed using the R package edgeR.The averages of three biological replicates (log 2 converted) are shown in the figures.